Boolean network model predicts cell cycle sequence of fission yeast 
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A Boolean network model of the cell-cycle regulatory network of fission yeast (Schizosaccha- 
romyces Pombe) is constructed solely on the basis of the known biochemical interaction topology. 
Simulating the model in the computer, faithfully reproduces the known sequence of regulatory activ- 
ity patterns along the cell cycle of the living cell. Contrary to existing differential equation models, 
no parameters enter the model except the structure of the regulatory circuitry. The dynamical 
properties of the model indicate that the biological dynamical sequence is robustly implemented 
in the regulatory network, with the biological stationary state Gl corresponding to the dominant 
attractor in state space, and with the biological regulatory sequence being a strongly attractive 
trajectory. Comparing the fission yeast cell-cycle model to a similar model of the corresponding 
network in S. cerevisiae, a remarkable difference in circuitry, as well as dynamics is observed. While 
the latter operates in a strongly damped mode, driven by external excitation, the S. pombe network 
, represents an auto-excited system with external damping. 
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tions; robustness 

Introduction 

Predicting the dynamics of complex molecular networks that control living organisms is a central challenge of 
systems biology. While cell-wide, or organism- wide, models of genetic and molecular interactions appear well out of 
reach, predictive models of single pathways and small modular molecular networks of living cells have been studied 
with great success and are a matter of active research [J [H, H, El ■ 

Given that the biochemical details of a chemical molecular network are known, standard techniques are at hand for 
their computer simulation. A method capturing molecular details is to use chemical Monte-Carlo simulations [H, [||, 
less computationally costly and perhaps the most commonly used approach to modeling biochemical pathways and 
networks are differential equations which capture the underlying reaction kinetics in terms of rates and concentrations 
Q. This method is highly developed today and is broadly applied to predictive dynamical modeling from single 
pathways to complex biochemical networks [8|. 

Such mathematical models contain detailed information about the time evolution of the system which, in some 
circumstances, is more than we are interested in. For many biological questions, knowledge of the sequential pattern 
^vq ■ of states of the central control circuit of a cell would be a sufficient answer, as, for example, in cell cycle progression, 
cell commitment (e.g. to apoptosis), and in stem cell control and differentiation. When we are interested in the path 
that a cell takes, the exact time course of the control circuit dynamics may not be needed, however, its modeling 
' takes most effort and often one needs to know large numbers of biochemical parameters that are not easily obtained 

Indeed, recent research indicates that some molecular control networks are so robustly designed that timing is not 
a critical factor Vice versa, as a working hypothesis, this observation bears the chance for vastly simplified 

■ dynamical models for molecular networks, as soon as one drops the requirement for accurate reproduction of timing 
by the model, just asking for the sequence of dynamical patterns of the network. Recent studies demonstrate, that 
such more simplified models indeed can reproduce the sequence of states in biological systems. For example, a class 
of discrete dynamical systems with binary states, mathematically similar to models used in artificial neural networks, 
has recently proven to predict specific sequence patterns of expressed genes as observed in living cells [HI, [l3| . 

Such models are in the mathematical tradition of random Boolean networks which, for decades, served as a simplistic 
analogy for how gene regulation networks could in principle work [l4|- I n these historical studies, dynamical properties 
of random networks of discrete dynamical elements were studied to derive possible properties of (the then hardly 
known) regulatory circuits [l5j . In the new approach outlined above, however, similar mathematical elements now 
serve to simulate one specific known biological control network. From a different perspective, they can be viewed as 
a further simplification of the differential equation approach [l6| . Recent application of this model class to modeling 
real biological genetic circuits show that they can predict expression pattern sequences with much less input (e.g. 
parameters) to the model as the classical differential equations approach. Examples are models of the genetic network 
of A. thaliana [l7l. llH. [l9| . the cell-cycle networks of S. cerevisiae [HI and of the mammalian cell-cycle [2(j, as well as 
the segment polarity gene network in D. melanogaster [l2l . l2l| . 

For example, the model by Albert and Othmer [12[ of the segment polarity gene network in D. melanogaster, as well 
as the model by Li et al. [13l | of the S. cerevisiae cell-cycle control network, yield accurate predictions of sequential 
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expression patterns, previously not obtained from such a simple model class. In these models, the dynamics can be 
viewed in terms of flow in state space of possible states of the network, converging towards so-called attractors, or 
fixed points, which here correspond to specific biological states. These attractors and their basins of attraction in 
state space mainly depend on the circuitry of the network, and their analysis yields further information about the 
robustness of the dynamics against errors or mutations. 

How generic is this approach? In this article we address the question whether the approach of discrete dynamical 
network models is a more general method, namely whether constructing predictive dynamical models for gene regula- 
tion from Boolean networks is a straightforward procedure that generalizes to other organisms. We choose the fission 
yeast (Schizosaccharomyces Pombe) cell-cycle as an example system that on the one hand is well understood in terms 
of conventional differential equation models, but on the other hand is markedly different from the above examples, as 
S. cerevisiae. S. Pombe has been sequenced in 1999 and has been used as a model organism only relatively recently 
[22| . Models exist [H, [24[ that mathematically model the fission yeast cell-cycle with a common ODE (ordinary 
differential equation) approach. These are based on a set of differential equations for the biochemical concentrations 
that take part in the network and their change in time (and space) . This approach allows to predict the dynamics of 
the fission yeast cell-cycle for the wild- type and some known mutant cells [13, [2{| . 

We will in the following construct a discrete dynamical model for the fission yeast cell cycle network. An interesting 
question will be, how far we will get without considering parameters, as kinetic constants etc., that are a key ingredient 
of the existing models. We will base our model on the circuitry of the known biochemical network, only. Let us in 
the next section briefly review the fission yeast cell cycle network, then define our discrete dynamical model in the 
subsequent section. This is followed by a section reporting our results, and then we will compare our findings with a 
similar model of the budding yeast (S. cerevisiae) network and conclude with a discussion. 



Let us briefly review the regulatory processes that control the cell cycle in Saccharomyces Pombe. The full process 
of one cell division consists of four stages, named G1-S-G2-M. At the first stage (Gl), the cell grows and, under 
specific conditions, commits to division. At the second stage (S), DNA is synthesized and chromosomes are replicated. 
This is followed by a "gap" stage G2. The final stage (M) corresponds to mitosis, in which chromosomes are separated 
and the cell divides itself. Eventually, after the M stage, the cell enters Gl again, thereby completing one cycle. 

The biochemical reactions that form the network that controls the fission yeast cell-cycle have been studied in 
detail over the last years [2(| H3, [H, [29|, [3(| Hi], HH H3 • The major role is played by a cyclin-dependent protein 
kinase complex Cdc2/Cdcl3 with a protein Tyr-15, a residue of Cdc2. Tyr-15 acts as a label for high Cdc2/Cdcl3 
concentration. It is inactive during the G2 phase, when Cdc2/Cdcl3 is phosphorylated, and becomes active during 
the G2-M transition 0, [25| . The other members that participate in the cell-cycle control can be attributed to two 
different classes. The first class consists of positive regulators of the kinase Cdc2/Cdcl3: "Start kinase" (SK), a group 
of Cdk/cyclin complexes (Cdc2 with Cigl, Cig2 and Pucl cyclins), and the phosphatase Cdc25. A second class is 
composed of the antagonists of the complex Cdc2/Cdcl3: Slpl, Ruml, Stc9, and the phosphatase PP 0. 

We give a full compilation of the network of key-regulators of the fission yeast cell cycle network in Table 1, 
corresponding to our current knowledge as given in 

nnnm Ais ° 

our translation into an interaction graph with 
activating and inhibiting links is given in the table, which is the starting point for our discrete dynamical network 
simulation of this network. Let us in the next section define the discrete dynamics that we will simulate on this graph. 



We assume proteins to be the nodes of the network and assign a binary value Si(t) £ {0, 1} to each node i, denoting 
whether the protein is present or not (due to different possible biochemical mechanisms, as, e.g., gene expression of 
a corresponding protein, or fast biochemical reactions as phosphorylization). The interactions between the nodes, as 
compiled in Table 1, are denoted as links, or arrows, see Figure 1. We do not quantify any interaction strength, except 
whether a link is present or not, and whether it is activating or inhibiting. Again, different biochemical mechanisms 
are subsumed under this simplified picture, as, e.g., transcriptional regulation, or faster enzymatic interactions. The 
dynamics of the nodes are updated (in parallel) in discrete time steps according to the following rule: 



The fission yeast cell cycle network 



A discrete dynamical model of the cell cycle network 
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TABLE I: The rules of interaction of the main elements involved in the fission yeast cell cycle regulation. 



where a,ij = 1 for an activating interaction (green link) from node j to node i, and a,j = — 1 for an inhibiting (red) 
link from node j to node i, and a*j = for no interaction at all. This definition follows closely the approach in [l3|. 
The dramatic simplification steps in constructing this model consist in not differentiating between absolute values of 
interaction strengths on the one hand, and not distinguishing between the different time scales of the biochemical 
interactions involved on the other. This corresponds to dropping all biochemical parameter values, time constants as 
well as binding constants, from the differential equation models. As we will see below, dynamical models on networks 
can be built to be insensitive to these parameters, provided that the interaction topology has certain properties. 

Two of the ten proteins included in the model exhibit a slightly different activation behavior, which we account 
for by the two following rules (which alternatively could be incorporated into the above equations as a non-zero 
activation threshold). Slpl is only activated by a highly active complex Cdc2/Cdcl3, which corresponds to both 
active Cdc2/Cdcl3 as well as active Tyrl5, since Tyrl5 labels the level of activity of Cdc2. This mechanism acts as 
a barrier for entering mitosis. The second special rule is to add "self-activation" (corresponding to adding a negative 
activation threshold) to the node Cdc2/Cdcl3, as it is otherwise not positively regulated. 

We also follow [l3[ by adding " self-degradation" (yellow loops) to those nodes that are not negatively regulated by 
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FIG. 1: Network model of the fission yeast cell-cycle regulation. 



others, representing the continuous degradation of proteins in the cell, which corresponds to an = — 1. 

Nodes, that have the same function as, for example, Weel/Mikl and SK (Cdc2/Cigl, Cdc2/Cig2, Cdc2/Pucl) are 
joined together in a single node (see Figure 1), as it does not make a difference in the specific mathematical model 
dynamics considered here. 

Finally let us define the initial condition of the model at the start of the simulation, which is chosen to correspond 
to the biological start condition, i.e. all nodes being in the OFF (inactive) state, except for the proteins Start, Ste9, 
Ruml, and Weel/Mikl [HI]. 

Results 

Simulation of the fission yeast cell cycle 

Let us first consider the time evolution of the proteins of the dynamical model described above. Let us run the 
cell-cycle model by exciting the Gl stationary state with the cell size signal ("Start" node). This initiates a sequence 
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TABLE II: Temporal evolution of protein states in the cell cycle network. 



of network activation states ("expression patterns") that, eventually, return to the Gl stationary state. The temporal 
evolution of the protein states is presented in Table 2, where one observes a sequence of states which exactly matches 
the corresponding biological expression pattern along the cell-cycle, from the excited Gl state (START) through S 
and G2 to the M phase and finally back to the stationary Gl state. It corresponds to the biological time sequence 
of the protein states in the cell-cycle control network. This is a remarkable observation as it is unlikely to occur by 
chance due to the size of the state space. 

In the next step we run the model starting from each one of the 2 10 = 1024 possible initial states. We find that 
each initial state flows into one of 15 stationary states (fixed points). The largest attractor belongs to a fixed point 
attracting 77% of all network states. Our first observation is that this fixed point exactly coincides with the biological 
Gl stationary state (see Table 3) of the cell. Thus, the biological target state is the dominant attractor of the network 
dynamics. As soon as the system reaches this state with the specific corresponding combination of active and inactive 
proteins, it stays there, and is likely to do so even in the presence of perturbations. 

A further observation is best depicted by Figure 2, showing the dynamical flow of the network states, and how it 
converges towards the biological fixed point. In this figure, the dynamical trajectories in the state space starting from 
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TABLE III: All attractors (fixed points) of the dynamics of the network model for the fission yeast cell cycle regulation. 




FIG. 2: State space of the 1024 possible network states (green circles) and their dynamical trajectories, all converging towards 
fixed point attractors. Each circle corresponds to one specific network state with each of the ten proteins being in one specific 
activation state (active/inactive). The largest attractor tree corresponds to all network states flowing to the Gl fixed point 
(blue node). Arrows between the network states indicate the direction of the dynamical flow from one network state to its 
subsequent state. The fission yeast cell-cycle sequence is shown with blue arrows. 
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TABLE IV: Homologue proteins related to the cell cycle networks of fission yeast and budding yeast 



all 1024 possible initial states of the network are shown. Each network state is represented by a dot, with the arrows 
between them indicating the dynamical transition from one state to its temporally subsequent state. At the root of 
the largest attractor (tree) the Gl state is found and the blue arrows show the biological time sequence that leads to 
it. This attractor tree consists of 77% of all network states. 

We further performed a robustness test by reversing the state of a single, randomly chosen node while the network 
proceeds through the biological sequence. This deviation from the biological pathway by the activity state of one 
single protein at one randomly chosen step of the cycle, the system returns to the fixed point Gl in 90 out of 100 
possible cases. Thus we observe an additional robustness in the fission yeast cell-cycle network, meaning that there is 
an increased probability to stay in the attractor basin of the biological fixed point when perturbing states along the 
biological trajectory. 

An immediate question about the specific network structure considered here is whether the architecture of the 
network has special properties as, for example, traces of being optimized by biological evolution. We compare the 
network dynamics to the null model of random networks with the same number of inhibiting and activating links, 
self-degrading and self- activating nodes and the same activation thresholds. Indeed one finds that the corresponding 
random networks typically have smaller attractors. The mean size of the biggest attractors is about 38% of all initial 
states (averaged over 1000 random networks). This may indicate that attractor basin size of the biological attractor 
is optimized, possibly in order to provide additional dynamical robustness. 

Comparison with S. cerevisiae 

The two yeasts, S. cerevisiae and S. pombe, are remarkably different cells and a comparison may provide insights 
relevant for the understanding of higher eukaryotic organisms. As we now have discrete dynamical models for the cell 
cycle network of both of them at hand (this work, as well as [HI), let us discuss how they compare. 

As these two organisms are closely related genetically, one might expect a large overlap also in the biochemical 
control machinery. On the other hand, the biology of the two is markedly different, so there have to be some differences 
on the biochemical level as well. As an overview, the second model is shown in Figure 3. 

There are a number of closely related genes (see Table 4) between the two yeasts [lj| which, however, can have 
vastly differing functions [22J. In fission yeast, for example, phosphatase Cdc25 is required for the G2-M transition, 
while in the model of budding yeast [l3[ the corresponding homologue Mihl is insignificant. The reason is that in 
the fission yeast cell cycle, Cdc25 removes an inhibitory phosphate group from the residue Tyr-15 of Cdc2, which 
is important for the right timing of the G2-M transition. In contrast, the tyrosine residue in S. cerevisiae Cdc28 
kinase (fission yeast: Cdc2) is not as critical and usually not phosphorylated. Therefore, for a model of fission yeast, 
Cdc25 is essential, whereas the homologue Mihl in budding yeast is not [l3| . One other example is the role of the 
protein Cdcl3. In fission yeast it acts in a complex with Cdc2, while in the budding yeast model its functionality 
is represented by two complexes Clbl,2/Cdc28 and Clb5,6/Cdc28, which exhibit some differences in interactions, as 
well as in timing. 

Despite of the differences in many details, the general logic of both yeast cell cycles is surprisingly similar and 
exhibits a number of " structural homologues" . For example both exhibit a negative feedback loop similar in role: 
Clbl,2/Cdc28 activates Cdc20 which inhibits Clbl,2/Cdc28 (fission yeast: Cdc2/Cdcl3, Tyrl5 activate Slpl,which 
inhibits Cdc2/Cdcl3). 

The most interesting comparison is in our view on the level of the global network dynamics. From this point of 
view, the S. cerevisiae network is a strongly damped system, driven by external excitation. External signals are 
entering the network, triggering signal cascades in the network that induce the subsequent phases. In contrast, the 
network of 5". pombe corresponds to an auto-excited system (there are two nodes with self-excitation - Cdc2/Cdcl3 
and Weel/Mikl) with additional damping. Here, an external signal works as a trigger mechanism that counteracts 
internal damping, causing the auto-excitation to spread its activity in the system. 

While these differences in the "mechanics" of the signalling networks are considerable, the overall dynamics is 
surprisingly similar. The state space picture is quite similar in both cases: one observes only a small number of 
attractors and just one big global attractor (with 86% resp. 77% of all initial states) which for both organisms 
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FIG. 3: Budding yeast cell cycle network model of [TjJ, for comparison with our model of fission yeast. This network relies 
more on transcriptional regulation than the fission yeast network (note that some homologues corresponding to the latter do 
not have to be included here). Note also the difference in circuitry. 

corresponds to the stationary Gl state. 

Finally, a most prominent difference between the two yeast networks is their choice in biochemical machinery: S. 
cerevisiae relies more on transcriptional factors while S. pombe mostly relies on post-translational regulation [34| . 
From the methodological point of view, we note that for this reason we were surprised to find our model for the S. 
pombe cell cycle network so robust against neglecting the vastly different time scales of interactions, which we expected 
to be the major difficulty in constructing a discrete dynamical model for S. pombe as compared to S. cerevisiae. 

Discussion 

We have constructed a Boolean model for the biochemical network that controls the cell cycle progression in fission 
yeast S. pombe, and found a number of interesting results. The dynamics of this network reproduces the time sequence 
of expression patterns along the biological cell cycle, solely on the basis of the connectivity graph of the network, 
neglecting all biochemical kinetic parameters. The dynamics of the network is characterized by a dominant attractor 
in the space of all possible states, with an attractor basin that attracts most of all states. The network dynamics are 
robust against perturbation of the biological expression pattern. 

The results obtained from our model are in accordance with the existing ODE model of fission yeast [l(|. Let us 
discuss the differences between these two approaches. The S. pombe ODE system [l(| has several steady state solutions. 
One can identify every such solution with the corresponding physiological stage. The growth of cell size brings the 
cell from one phase to another via a series of bifurcations. At the same time, other variables indicate the degree of 
activity of various components of the cell regulatory nodes. One observes [25[ that the typical curves depicting this 
activity have almost rectangular shape. This motivates our choice of binary valued function to approximate protein 
concentrations in time. Further, the ODE-bascd model makes use of continuous system parameters, which we omit 
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and replace by their signs, only. As a result, the ODE bifurcation curve then corresponds to the Boolean biological 
path. The main advantage of our Boolean model is that we were able to drop 47 kinetic constants that were necessary 
in the ODE approach and, while doing so, still reproduce the biological activation pattern of the system. 

This fact and our further observations point at built-in dynamical robustness of the network, which may provide 
a further mechanism for organisms to ensure functional robustness [35j . In return, our study indicates that the 
regulatory robustness of biological chemical networks may allow for "robust" modeling approaches: Our paradigm 
here is nothing but assuming that biochemical networks are functioning in a parameter-insensitive way — which 
motvated us to eliminate all tunable parameters from the model. That our model reproduces the biological sequence 
instantly without any further parameter tuning, confirms our assumption a posteriori. We therefore encourage further 
modeling experiments with the here presented quite minimalistic approach, as it may prove a quick approach to 
predicting biologically relevant dynamical features of genetic and protein networks in the living cell. 
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